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A STUDY OF THE MAIN RESONANCES OUTSIDE THE 
GEOSTATIONARY RING 


ALESSANDRA CELLETTI AND CATALIN GALES 

Abstract. We investigate the dynamics of satellites and space debris in external res¬ 
onances, namely in the region outside the geostationary ring. Precisely, we focus on the 
1:2, 1:3, 2:3 resonances, which are located at about 66 931.4 km, 87 705.0 km, 55 250.7 
km, respectively. Some of these resonances have been already exploited in space mis¬ 
sions, like XMM-Newton and Integral. 

Our study is mainly based on a Hamiltonian approach, which allows us to get fast 
and reliable information on the dynamics in the resonant regions. Significative results 
are obtained even by considering just the effect of the geopotential in the Hamilton¬ 
ian formulation. For objects (typically space debris) with high area-to-mass ratio the 
Hamiltonian includes also the effect of the solar radiation pressure. In addition, we 
perform a comparison with the numerical integration in Cartesian variables, including 
the geopotential, the gravitational attraction of Sun and Moon, and the solar radiation 
pressure. 

We implement some simple mathematical tools that allows us to get information on 
the terms which are dominant in the Fourier series expansion of the Hamiltonian around 
a given resonance, on the amplitude of the resonant islands and on the location of the 
equilibrium points. We also compute the Fast Lyapunov Indicators, which provide a 
cartography of the resonant regions, yielding the main dynamical features associated 
to the external resonances. We apply these techniques to analyze the 1:2, 1:3, 2:3 
resonances; we consider also the case of objects with large area-to-mass ratio and we 
provide an application to the case studies given by XMM-Newton and Integral. 


1. Introduction 

Resonances play a major role in the dynamics of satellites around the Earth; indeed, 
most of the satellites in MEO and GEO regions * 1 are positioned in correspondence with 
the 2:1 and 1:1 gravitational resonance (see, e.g., [6], [10], [12], [19], [20], [24]). Being 
in such resonances implies that the satellite makes two orbits or one orbit during one 
rotation of the Earth around its spin-axis. Other resonances might be important as 
well, and the aim of this work is to investigate some external resonances, precisely the 
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1:2, 1:3, 2:3 resonances, whose distances from the center of the Earth are, respectively, 
about 66 931.4 km, 87 705.0 km, 55 250.7 km. A remarkable fact is that space missions 
have already used these resonances; in particular, Integral (International Gamma-Ray 
Astrophysics Laboratory) has a semimajor axis corresponding to the 1:3 resonance, while 
the semimajor axis of XMM-Newton (X-ray Multi-Mirror Mission) is at the 1:2 resonance. 
We claim that such resonances can be particularly useful when dealing with space debris; 
indeed, the resonant dynamics can be exploited to move space debris in safe regions, either 
placing them in the stable equilibria, which prevent chaotic variations of semimajor axis, 
or moving the debris along the chaotic invariant manifold associated to the hyperbolic 
equilibria. 

In order to perform such investigation, we make use of the Hamiltonian formalism, 
allowing us to make a fast and accurate description of the dynamics. In particular, 
we introduce a model including a suitable expansion of the geopotential in spherical 
harmonics. The expansion is limited to a small number of terms, which are chosen by 
carefully evaluating which of them are the most significative ones in specific regions of 
the orbital parameter space. In this way, we reduce very much the computational time, 
although the model still retains the essential features of the dynamics. For objects with 
high area-to-mass ratio we include also a suitable expansion of the contribution due to 
the solar radiation pressure. Also in this case we use the Hamiltonian approach. 

This strategy allows us to obtain several information on the dynamics within a very 
short computational time, most notably the location of the equilibrium points, their 
dependence on the orbital parameters, the size of the resonant regions. Further informa¬ 
tion are obtained by making a cartography ([8]) of the different regions by means of the 
computation of the Fast Lyapunov Indicators ([7], see also [1]). In this way we obtain 
very detailed maps, showing the long term evolution of the semi-major axis, although 
within the Hamiltonian approach we limit the cartography to a model including just the 
geopotential and, possibly, the solar radiation pressure. As a consequence, to validate 
our results we make a comparison with the Cartesian equations of motion by computing 
maps which include also the effects of Sun and Moon. Attention must be paid in com¬ 
paring the results and, precisely, while integrating the Cartesian equations of motion, we 
need to transform from osculating to mean orbital elements. 

We also provide an application of our technique to two sample cases: XMM-Newton, 
which is related to the 1:2 resonance, and Integral, related to the 1:3 resonance. These 
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space missions are characterized by a large eccentricity; henceforth, a dedicated expansion 
of the geopotential is necessary, in order to include terms which are relevant at high 
eccentricities. Again, we use the Cartesian approach to validate the results obtained 
through the Hamiltonian formalism, although we find that in the case of XMM-Newton 
the effect of the Moon shapes the resonant islands and should be included in the overall 
discussion to obtain an accurate description of the dynamics. 

This paper is organized as follows. In Section 2 we introduce the Cartesian equations 
of motion, including the geopotential, the influence of Sun and Moon, and the solar 
radiation pressure. In Section 3 we introduce the Hamiltonian of the geopotential and we 
provide explicit expressions for the secular and resonant expansions of the Hamiltonian. 
In Section 4 we make a qualitative analysis of the resonant regions, by reducing the 
study to a very limited number of terms of the expansion and by using the pendulum¬ 
like structure of the resonant zone to compute the amplitude of the resonant islands. 
A cartography based on the computation of the Fast Lyapunov Indicators is given in 
Section 5, while the case of large area-to-mass ratio is investigated in Section 6. The 
cases of XMM-Newton and Integral are studied in Section 7, where the expansions of 
the geopotential have been extended to encompass the case of large eccentricities. Some 
conclusions and perspectives are presented in Section 8. 


2. Cartesian and Hamiltonian equations of motion 


We consider a small body, say S, moving in the gravitational held of the Earth, which 
we assume to be oblate, and subject to the influence of Sun, Moon and to the effect of 
the solar radiation pressure (hereafter SRP). We assume that the body is so small, that 
its influence on Earth, Sun and Moon can be neglected. Let r = (x, y, z) be the radius 
vector of S in a geocentric quasi-inertial fixed frame. 

The corresponding equations of motion are the sum of the equations describing the 
Earth’s gravitational influence, the oblateness effect, the solar and lunar attraction, and 
the SRP. Let mg, him be the masses of Sun and Moon, rg, r M the position vectors of 
Sun and Moon (whose explicit expressions are given, e.g., in [18]), Q the gravitational 
constant. The equations of motion of S can be written as 
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where R 3 is the rotation matrix about the Earth’s polar axis, 6 is the sidereal time, V is 
the gradient in the synodic frame, V (r) is the force function due to the attraction of the 
Earth: 


E(r) =g f 

JVe 


pOp) 


(JVe 


where p(r p ) is the density at some point r p in the Earth and Ve is the volume of the Earth. 
The second and third terms in (2.1) model the attraction of Sun and Moon, respectively. 
The last term in (2.1) models the SRP and depends on the adimensional reflectivity 
coefficient C r (fixed to 1 in this paper), the radiation pressure P r = 4.56 • 10 -6 N/m 2 for 
a body located at as = 1 AU, the area-to-mass ratio A/m with A the cross-section of 
S and m its mass. 


We proceed now to write in Hamiltonian formulation the term in (2.1) corresponding 
to the geopotential. To this end, we introduce the action-angle Delaunay variables, 
denoted as (L, G, H, M, to, fl), where L = ^//i B a, G = L\J\ — e 2 , H = Gcosi with a 
the semimajor axis, e the eccentricity, i the inclination, M the mean anomaly, oo the 
argument of perigee, D the longitude of the ascending node and Pe = Gm E with the 
mass of the Earth. The geopotential Hamiltonian can then be written as (see [2]) 

U{L , G, H, M, u, n, 0) = + R ear th(L , G, H, M, cn, 0, 9) , (2.2) 

where R ea rth represents the perturbing function, whose explicit expression is given as 
follows. 

In the geocentric quasi-inertial frame, the geopotential in (2.2) can be expanded as 

([ 11 ]) 

00 n D n 00 

Rearth = £>nmp(i) ^ S nmpq (M, CO, Q, 9) , (2.3) 

n =2 m =0 p =0 q=—oo 

where the inclination and eccentricity functions F nmp , G npq can be computed by well- 
known recursive formulae (see, e.g., [11]). The angle S nmpq is defined as 
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where C nm and S nm are, respectively, the cosine and sine coefficients of the spherical 
harmonics potential terms (see Table 1 for concrete values) and 


4/ nmpq = (n — 2 p)co + (n - 2 p + q)M + m(fl - 6 ) . 


(2.5) 
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As common in geodesy, we introduce also the quantities J nm defined by 
Jnm = \J A nm K TTl 0 , J n Q = J n = CnO 
and the quantities X nm defined through the relations 

Cnm Jnm COS(777. Xnm) ; Rnm Jnm sin(/7lA rlm ) . 


3. Secular and resonant Hamiltonian 

The expansion of the function R ear th in (2.3) contains an infinite number of trigonomet¬ 
ric terms, but the long term variation of the dynamics is mainly governed by the secular 
and resonant terms. Let us recall that a gravitational resonance (also called a tesseral 
resonance) occurs whenever there is a commensurability between the orbital period of 
the minor body and the period of rotation of the Earth. 

In astronomical applications the condition of gravitational resonance is satisfied only 
within a certain approximation and cannot be obviously satisfied exactly. 

Notice that, by using Kepler’s third law, a j : i resonance corresponds to a semimajor 
axis equal to a^x = a geo , where a geo = 42164.1696 km is the semimajor axis 

corresponding to the geostationary orbit. From this formula we derive the location of 
the 1:2 resonance at 66 931.4 km, the 1:3 resonance is found at 87 705.0 km, while the 
2:3 resonance is at 55 250.7 km. 


Let us expand the Earth’s gravitational potential up to terms of degree and order 
n = m = N for some iV > 0; then, we approximate R ea rth with the expression 

N n n oo 

td _ jDsec rxres . -nnonres \ ^ \ ^ \ ^ \ ^ 

-dearth earth ' ^earth ' earth / v / j / j / j 'nmpq ? 

71=2 m= 0 p= 0 q=— oo 


where R^rthi J^ewrthi J^earth S denote, respectively, the secular, resonant and non-resonant 
contributions to the Earth’s potential and where we have introduced the coefficients: 


Jnmpq 


HeR 
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n+1 " nm P 


npqi^fJ) S n mpq(XJ\ UJ , 12 , 0 ) 


(3.1) 


Remark 1. From (2.5) we see that the quantity T nmpq depends also on 6 j, Fl, which can 
be small, but not exactly zero. As a consequence, the stationary solutions have different 
locations according to the values of the integers n, m, p, q. As already noticed in [2], this 
implies that each resonance might split into a multiplet of resonances. 
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Our next task will consist in the determination of the secular part of the expansion 
(2.3) by computing the average over the fast angles, say R^thi and the determination of 
the resonant part associated to a given j : l gravitational resonance, say R^lth- 

Since the value of the oblateness coefficient J 2 = J20 is much larger than the value 
of any other zonal coefficient (see Table 1), we consider the same secular part for all 
resonances; the explicit expression will be given in Section 3.1 and Appendix A. 


n 

m 

c 

w nm 

Q 

u nm 

Jnm 

^ nm 

2 

0 

-1082.6261 

0 

1082.6261 

0 

2 

1 

-0.000267 

0.0017873 

0.001807 

—81°5116 

2 

2 

1.57462 

-0.90387 

1.81559 

75°0715 

3 

0 

2.53241 

0 

-2.53241 

0 

3 

1 

2.19315 

0.268087 

2.20947 

186°9692 

3 

2 

0.30904 

-0.211431 

0.37445 

72°8111 

4 

0 

1.6199 

0 

-1.619331 

0 

4 

1 

-0.50864 

-0.449265 

0.67864 

41°4529 

4 

2 

0.078374 

0.148135 

0.16759 

121°0589 

4 

3 

0.059215 

-0.012009 

0.060421 

56°1784 


Table 1. The coefficients C nm , S nm , J nm (in units of 10 6 ) up to degree 
and order 5; values computed from [5] (see also [3], [18]). 


The expansions of the resonant part, say RfffJh, are composed by several different 
terms, say 7K for some k G Z + ; in practical computations it is essential to retain the 
minimum number of significant terms. To this end, we introduce the following heuristic 
definition of dominant term. 

Definition 2. Consider a j : t gravitational resonance and let \ <jf) be the associated 
stroboscopic mean node (see [9], [14],). Given the orbital elements ( a,e,i), we say that a 
term for some k G Z + of the expansion of Rfff r j^', say %k = gk{a, e, i) cos(k + 
7 k) for some 7 k constant, is dominant with respect to the other harmonic terms of the 
resonant part, if the size of \gk(a,e,i)\ is bigger than the size of any other term of the 
expansion. 

Making use of Definition 2, for a given j : i gravitational resonance, we proceed to 
approximate the Hamiltonian function by the expression 

njresj'.t _ - 1 jDsec jyres j\i 

O r 2 ' earth ' ^earth ’ 










A STUDY OF THE MAIN RESONANCES OUTSIDE THE GEOSTATIONARY RING 7 

where for the resonant part R^rtii we considered an optimal degree N of the expansion, 
obtained using the following algorithm: 

i) expand the resonant part up to a degree n large enough (at least n = N + 1); 

ii) plot the dominant terms of the expansion as a function of e, i within prescribed 
intervals that we choose to be e 6 [0, 0.5], i e [0, 90°]; 

Hi) compute the size of each function \gk(aj-.e, e, z)|; 

iv) through a color scale represent the index k of the dominant term in the plane e — i, 

v ) the optimal degree N of the expansion is the largest degree of the terms, whose 
index is represented in the plot of dominant terms computed in iv). 

For the 1:2 and 1:3 resonances the optimal degree of expansion of R^rth is N = 4, 
while for the 2:3 resonance the optimal degree is 3. In the next sections we provide the 
secular and resonant terms up to the second order in the eccentricity; explicit expressions 
of all coefficients providing the secular and resonant terms are given in Appendix A. 

3.1. The secular part of R ear th- With reference to the expression for S nmpq given in 
(2.4)-(2.5), the secular terms correspond to m = 0 and n — 2p + q = 0. We consider the 
secular part up to terms of degree n = 4 and neglect higher order harmonic terms. 

Since G npq (e) = O(e^), then up to second order in the eccentricity, the secular part is 
given by 

Rearth ~ 7200-2 + 72010 + 72022 + 7 - 301-1 + 73021 + 7401-2 + 74020 + 74032 , ( 3 - 2 ) 

where the functions T nm pq have been introduced in (3.1). Using the formulae given by 
[11] for the functions F nmp and G npq , the explicit expression of the above terms can be 
found in Appendix A. 


3.2. The resonant part of R ear th • From (2.4)-(2.5) we see that the terms associated 
to a resonance of order j : i correspond to j (n — 2p + q) = £m. We consider the resonant 
part up to degree and order n = m = N with IV = 4 for the 1:2, 1:3 resonances and 
IV = 3 for the 2:3 resonance. The terms whose sum provides the resonant part of R^rth 
for the 1:2, 1:3, 2:3 resonances are listed in Table 2, while their explicit expressions are 
given in Appendix A. 


Remark 3. Except for the 1:3 resonance, the expansions up to the second order in the 
eccentricity of R™*ih describe with a good enough accuracy the main dynamical features 
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j:£ 

N 

terms 

1:2 

4 

72100; 72112; 72202 ; 7310-1, 73111, 73201; 7410-2; 74110; 74122 ; 74200 ; 74212; 74302 

1:3 

4 

72101; 72204 ; 73100; 73112, 7410-1, 74111, 74202 

2:3 

3 

72201, 73200 , 73212 


Table 2. Terms whose sum provides RfffJth't ^ ie explicit expressions are 
given in Appendix A. 


of the resonances when the eccentricity is in the range e E [0,0.5]. However, for the 1:3 
resonance, there is a term of order four in the eccentricity, precisely T 2204 , which plays a 
very important role for moderate and large eccentricities (see Section 5.3). 

4. Qualitative analysis 

Following [2], making use of elementary mathematical methods, we perform a first order 
qualitative analysis to select the dominant terms in specific regions of the parameter space 
(Section 4.1). After this selection, we proceed to give an estimate of the amplitude of the 
resonant islands (Section 4.2). This strategy will allow us to obtain useful information 
on the dynamics of the resonances within a very limited computational time. 

4.1. Dominant terms. Starting from the series expansions provided in Sections 3.1 
and 3.2, for each value of e, i, we compute the dominant terms according to Definition 2. 
Hence, we obtain the plots for the 1:2, 1:3, 2:3 resonances, as shown in Figure 1, upper 
panels. The orbital elements (e,i) are taken within the intervals e G [0,0.5], i E [0°, 90°]. 

Since the magnitude of the resonant terms decays with the degree n as (i?g/a) n , and 
since the semimajor axis is very large for the considered resonances, the most relevant 
role is played by the harmonic terms of order (7(722), even though they are of high order 
in the eccentricity. From Table 3 we see that we are led to consider at most 5 terms, thus 
reducing considerably the expansions for a specific resonance. 


Resonance 

Black 

Yellow 

Brown 

Blue 

White 

1:2 

12202 

74110 

73201 

73111 

74200 

1:3 

72204 

73100 

73112 

- 

- 

2:3 

72201 

73200 

- 

- 

- 


TABLE 3. The colors used in Figure 1 to represent the dominant terms. 
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Resonance 1:2 



0 0.01 0.02 0.03 0.04 0.05 

e 



0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 

e 



Resonance 1:2 Resonance 1:3 



Resonance 2:3 



0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 


e 


FIGURE 1. Left 1:2, center 1:3, right 2:3 resonance. Upper plots: 
dominant terms; for the 1:2 and 2:3 resonances we magnified the region 
[0, 0.05] x [0°, 90°], because in the rest of the domain the color is black. 
The colors used to represent the dominant terms are explained in Table 3. 
Lower plots: amplitude of the resonances for co — 0°, Q = 0°; the color 
bar provides the measure of the amplitude in kilometers. 


4.2. Amplitude of the resonant islands. The amplitude of the island around a given 
j : i resonance is obtained by making an expansion of the Hamiltonian around the 
resonance, thus reducing the problem to a pendulum-like Hamiltonian. This expansion 
will provide an analytical estimate of the amplitude of the resonant region (we refer to 
[2] for full details). Let us briefly summarize the method devised in [2], We start by 
taking into account the secular part and, precisely, just the largest term of the resonant 
part. Let the resonant Hamiltonian be written as 


n resj: \L, G, H, £M—j9, co, fi) = 


hi 


2 L 2 


±R S e e artk(L, G, H, u)+I%%*(L, G, H , £M-j9, co, fi) 

(4.1) 


with 


RZtt(L,G,H,lM-je,u,n) = 


Ni N 2 N 3 


EEE R^\l, G, H)cs(ki(£M — j6) + k 2 co + k^fl) 


ki = l /c2 — 1 k$ = l 


for some integers Ni, N 2 , N 3 and for some Fourier coefficients R^’ e \ In the above relation, 
cs could be either cosine or sine and k = (fc 1; k 2 , k^). Let L = L res be the resonant value of 
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the action L , which can be evaluated through Kepler’s third law. Expanding (4.1) around 
L res up to second order and retaining only the largest term in the resonant Hamiltonian, 
we obtain a Hamiltonian of the form 


'Hmax' e (L, G, H, £M — j6, oj, O) = a(L - L res ) - /3(L - L res f 

+ 7 cs(k™ ax (m - jo) + k™ ax uj + k™ ax n ), 

for suitable functions which can be approximated as a ~ /i| ; /L^ es , (3 ~ 3/d 2 E /(2Lj es ), while 
7 is the largest coefficient, say 7 = ( L res , G, H) with index /c = ( kY 1mx , k^ ax , k^ ax ). 

One can show (see [ 2 ]) that the amplitude of the j : £ resonant island is given by 

5 ' 

P 

We report in Figure 1 (lower panels) the amplitudes of the 1:2, 1:3, 2:3 resonances as a 
function of the eccentricity (between 0 and 0.5) and the inclination (between 0° and 90°); 
in the computations we fixed oj = 0° and H = 0°. The color bar provides the size of the 
amplitude in kilometers. A comparison with the results of Section 5 shows a remarkable 
agreement between these figures and the results we will obtain plotting the FLIs (using 
a much longer computational time). 


A 2 /2y 

2Aa — — ( — + 2 L res 
He \ p 


4.3. Casting the results. In this section we collect the information based on the fol¬ 
lowing items, which have been investigated in the previous sections: 

(1) the dominant terms with reference to Figure 1, which provide the terms of the 
expansions in Sections 3.1 and 3.2 leading the dynamics (see Section 4.1); 

( 2 ) the amplitude of the resonances associated to each dominant term (see Sec¬ 
tion 4.2); 

(3) the location of the equilibrium points associated to the different dominant terms, 
according to the expansion in Section 4.2 (see also Appendix A). 

This analysis has several advantages: (?) it is based on analytical arguments, (ii) it 
does not need to integrate the equations of motion, (in) it can be carried out within a 
very limited computer time. Indeed, the following information will provide a fast tool to 
obtain a qualitative behavior, which will be confirmed and refined by the cartographic 
study performed in Section 5. 


1:2 resonance. The dominant terms are 72202, 73 m, 73201, 74 no, 74200, but the most 
important one is T 2202 , which corresponds to the biggest amplitude of a few tenths of 
kilometers in the region in which e > 0 . 1 . Since the trigonometric argument of T 2202 


A STUDY OF THE MAIN RESONANCES OUTSIDE THE GEOSTATIONARY RING 


11 


is 2(a 12 — c o — A 22 ) (see Appendix A), one gets in the (< 712 , a) plane a single resonant 
island in an interval of length of 180° on the 0 T 2 axis, with the stable point located at 
< 7 12 = A 2 2 + oj + 180°fi:, k G Z. In general, the amplitude of the resonant island is bigger 
as the eccentricity is larger and the inclination is smaller. 

All these information are consistent with those shown in Figure 2; in fact, the equilibrium 
is located at about A 2 2 — 75° (see Table 1), while the amplitude is in agreement with the 
lower plots of Figure 1. 

1:3 resonance. The dominant terms are 72204, 73100 , 73112 ■ Although 72204 is of 
order four in the eccentricity, it is dominant for a large region of the (e, f)-plane. As 
a consequence, for moderate and large eccentricities the stable equilibrium points are 
located at cr 13 = A 22 + 2u + 180°k;, kGZ. On the contrary, for small eccentricities and 
nonzero inclinations, 73100 is dominant and therefore the elliptic equilibrium points are 
located at <r 13 = A 31 — 180° + 360°k, kgZ. 

Again, we find a very good agreement with the results of Section 5 (see, e.g., Figure 3). 

2:3 resonance. The dominant terms are 72201 , 73200 ; the term 73200 is dominant for 
small eccentricities (say e < 0.02). The term 72201 is dominant almost everywhere else. 
The width of the resonance is bigger as the eccentricity is larger. The stable equilibrium 
points are located at cr 23 = 2A 2 2 + uj + 360°ac, k£Z, when 72201 is dominant, while they 
are at cr 23 = 2A 3 2 + 90° + 360°/y kgZ, when 73200 is dominant (compare with Figure 4). 


5. Cartography 


The goal of this section is to describe some dynamical features of the external res¬ 
onances by means of the computation of the FLIs, which are suitable chaos indicators 
allowing us to distinguish between resonant, stable and chaotic motions. In Section 6 we 
will pay special attention to objects with large area-to-mass ratios. 

5.1. Fast Lyapunov Indicators. Making reference to [7], we define the FLI, which 
corresponds to the largest Lyapunov characteristic exponent at a fixed time, say t = T. 
Precisely, consider the n-dimensional differential system 


x = f(x) , 


where x e M n . We write the variational equations as 
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where x is an n-dimensional vector. For the initial conditions x(0) G M n , v(0) G R n and 
a time T > 0, the FLI is defined as 

FLI(x(0), v(0),T) = sup log ||v(t)|| . 

0 <t<T 

We will present different results based on the FLIs, computed in the plane of coordi¬ 
nates or in the parameter plane; the value of the FLI will be marked by a color scale, 
where darker colors represent a regular dynamics (either periodic or quasi-periodic), 
while lighter colors are associated to chaotic motions. 

5.2. Cartography of the 1:2 resonance. For the 1:2 resonance, we have twelve terms 
defining the resonant part FFJarth °f the geopotential (see Table 2). 

Except for eccentricities less than about 0.02 (see Figure 1, top left), the term 72202 
is dominant in the rest of the e — i domain. For moderate and large eccentricities, its 
magnitude is much larger that the size of any other term. As a consequence, it is natural 
to expect that pendulum like plots are obtained in the (er 12 ,a) plane (see Figure 2, top 
panels). The stable equilibrium point is located at cri 2 = A 22 + oj — 75° + u>. For large 
inclinations and small eccentricities the term 74200 is dominant and the equilibrium is 
located at <7 i 2 = (2A4 2 — 180°)/2 = 31° (modulus 180°). 

Figure 2, middle panels, plots the FLI values as a function of inclination and semimajor 
axis in order to evaluate the width of the resonance for each value of the inclination and to 
give a hint on the dynamics inside the resonance. For a specific inclination, the amplitude 
can be determined by measuring the distance between the two points on the separatrix 
obtained as the intersection of the vertical line corresponding to that specific inclination 
and the branches visible on the plot. For small eccentricities (Figure 2, middle left), the 
amplitude of the resonance is very small (at most two kilometers) in comparison with the 
amplitude for moderate and high eccentricities (compare with Figure 2, middle right). 

The bottom panels of Figure 2 show the FLI values as a function of eccentricity and 
semimajor axis for small and large inclinations, thus providing an estimate of the width 
of the resonance in terms of the eccentricity for a fixed value of the inclination. 

From these graphs, it follows that the amplitude of the resonance grows in magnitude 
with the increase of the eccentricity and it reduces in magnitude with the increase of the 
inclination (see Figure 2, bottom and middle panels). 

5.3. Cartography of the 1:3 resonance. A similar analysis can be performed for the 
1:3 resonance, where we have several terms defining R^rth ( see Table 2), among which 
the most important ones are 73100 and 72204 (see Figure 1, top center panel). 
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Figure 2. FLI map for the 1:2 resonance for uj = 0°, = 0°. Top 

panels: e = 0.1, i = 10° (left), e = 0.2, i = 10° (right). Middle panels: 
FLI in the (z, a) plane for e = 0.005, a 12 = 31° (left), e = 0.5, a 12 = 75° 
(right). Bottom panels: FLI in the (e.a) plane for ai 2 = 75° and i = 20° 
(left), i = 70° (right). 


Since 73100 is of order 0(1) (while the other terms are (7(e)), for small eccentricities 
some pendnlnm like plots are obtained, provided the inclination is not zero. The stable 
equilibrium point is located at cri 3 = A 3 i — 180° = 7° (see Figure 3, top left). 

On the other hand, since 72204 (which contains cos(2(ai 3 — 2 c j — A 22 ))) is dominant 
for large eccentricities, one gets two stable points, located at ai 3 = A 2 2 — 75° and 
(T i 3 = A 2 2 — 180° = —105° (see Figure 3, top right). Notice that, although the term 
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FIGURE 3. FLI for the 1:3 resonance for to — 0°, = 0°. Top panels: 

e = 0.005, i = 70° (left); e = 0.3, i = 25° (right). Middle panels: FLI in 
the (i,a) plane for e = 0.005, a 13 = 7° (left), e = 0.5, = 75° (right). 

Bottom panels: FLI in the (e,a) plane for cri 3 = 75° and i = 20° (left), 
i = 70° (right). 

72204 contains a fourth power of the eccentricity, the width of the resonance is larger in 
comparison with that associated to 73100 • 

The middle panels of Figure 3 show the FLI values as a function of inclination and 
semimajor axis, while the bottom panels refer to the (e, a) plane. 

5.4. Cartography of the 2:3 resonance. For the 2:3 resonance, we have three terms 
defining R^th ( see Table 2): 72201 , %200 and % 2 i 2 - The harmonic term 73200 is dominant 
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FIGURE 4. FLI for the 2:3 resonance for oj = 0°, = 0°. Top panels: 

e = 0.005, i = 10° (left); e = 0.005, i = 70° (middle) e = 0.3, i = 10° 
(right). Middle panels: FLI in the (i,a) plane for e = 0.005, a 23 = 236° 
(left), e = 0.5, <723 = 150° (right). Bottom panels: FLI in the (e, a) plane 
for a 23 = 150° and i = 20° (left), i = 70° (right). 


for small eccentricities and nonzero inclinations, while 72201 is dominant for the rest of the 
domain. Pendulum like plots are obtained for any value of eccentricity and inclination. 
For small eccentricities the stable equilibrium point is close to <j 2 3 = 2A 2 2 + ^ — 150°+ ca, 
when the inclination is small (see Figure 4 top left), while it is located in the vicinity 
of <723 — 2 A 32 + 90° = 236° when the inclination is large (see Figure 4 top middle). For 
moderate and high eccentricities, the stable point is located at a 2 3 = 2A 22 + oj = 150°+ ca 
(Figure 4 top right). 
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The amplitude of the resonant islands for small and large eccentricities can be deduced 
from Figure 4, middle panels, while the bottom plots provide the FLI as a function of 
eccentricity and semimajor axis for (relatively) small and large inclinations. 

5.5. A comparison with the full Cartesian equations. The purpose of this section 
is to complement the study performed by using the Hamiltonian formulation with some 
results obtained in Cartesian coordinates. In fact, the Hamiltonian formulation provides 
an approximation of the dynamics, due to the fact that we just consider the secular and 
resonant expansions. Moreover, these expansions have been developed up to second order 
in eccentricity. Clearly, the Cartesian equations of motion contain the non-truncated 
terms of the geopotential, including the short periodic terms, as well as the gravitational 
attraction of the Sun, Moon and the effect of the solar radiation pressure (see (6.1) 
below). 

Henceforth, the Cartesian formulation allows one to have a more complete model, 
although at the expense of a much higher computational time, at least with respect to 
the Hamiltonian formulation. However, the Cartesian approach cannot provide a global 
explanation of the resonant dynamics as a function of eccentricity and inclination, as it 
was done in the previous sections using the Hamiltonian formalism. The goal of this 
section will be to validate the results obtained by means of the resonant Hamiltonian 
describing an approximation of the geopotential (see Sections 5.2, 5.3, 5.4) by comparing 
such results with a more complete model, obtained using the Cartesian formalism and 
including all major effects. 

Let us remark that within the Hamiltonian formulation, we removed the short periodic 
perturbations by averaging over the fast angles, thus leading us to compute the mean 
orbital elements. As a consequence, when dealing with the equations of motion in Carte¬ 
sian coordinates, in order to represent the FLI as a function of the same variables as 
in the Hamiltonian approach, we need to transform from osculating orbital elements to 
mean elements. This computation implies a numerical average of the osculating elements, 
which is performed in the course of the integration itself; we refer to [25] for alternative 
approaches to compare osculating and mean orbital elements, based for example on nu¬ 
merical averaging, fixed-point iterations, or a Hamiltonian transformation. 

The results obtained by using the Hamiltonian formulation are validated by integrating 
the Cartesian equations of motion as shown in Figures 5. Moreover, we evaluate the 
effects of the other perturbing forces, in particular the influence of the Moon. We have 
used as starter a single step method (a Butcher numerical algorithm), while a multistep 
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predictor-corrector numerical method (Adams-Bashforth 12 steps and Adams-Moulton 
11 steps) performs most of the propagation. For each resonance, the total time span was 
20 000 sidereal days. For each orbit we used a fixed initial tangent vector; a different 
choice of the tangent vector might alter the results ([7], [15]), although it is not really 
relevant in our context and, indeed, we leave the exploration of the effect of a random 
choice of this vector to a future work. 

Comparing the above plots with those obtained using the Hamiltonian approach, we 
are led to the following conclusions: the main dynamical features of the 1:2, 2:3 and 1:3 
resonances, namely the location of the equilibrium points and the amplitude of resonant 
islands, which we already found by using the Hamiltonian formalism, are confirmed by 
integrating the full equations of motions. However, the other disturbing functions, in 
particular the perturbations due to the Moon, induce in some cases a substantial effect. 

In fact, from the plots of Figure 5 we notice that within the present range of eccen¬ 
tricities below 0.5 the influence of the Moon becomes larger for the orbits located farther 
from the Earth. For example, for the 1:2 resonance, only small variations are observed 
when adding the effects of Sun, Moon and SRP (see Figure 5, top left and middle panels) 
with a good agreement with the Hamiltonian formulation (compare with Figure 2, top 
right panel). Also the agreement with the 2:3 resonance is very satisfactory (compare 
Figure 5, top right with Figure 4, top right). On the contrary, in the case of the 1:3 
resonance, which is closer to the Moon, the Moon’s perturbation leads to large chaotic 
regions (see Figure 5, bottom panels and compare Figure 3, top right with Figure 5, bot¬ 
tom right). We will see in Section 7 that for large eccentricities the effect of the Moon is 
more important for the 1:2 resonance, than for the 1:3 resonance. The lunar interaction 
with the external resonances, as the orbital parameters (i.e., eccentricity and inclination) 
are varied, will be the subject of a further investigation. 

6. Dynamics of large area-to-mass ratio objects 

Extended analyses of the dynamics of high area-to-mass ratio objects are provided in 
several papers (see, e.g., [4], [16], [17], [21], [22], [23]). The purpose of this section is 
to explore the interaction of such objects with the 1:2, 2:3, 1:3 resonances, by using the 
Hamiltonian formalism. For simplicity, we focus our attention on planar orbits, for which 
we compute the FLI for various values of the area-to-mass parameter. To this end, we 
proceed to expand the disturbing function corresponding to SRP in terms of the orbital 
elements. Precisely, the potential describing the effect of SRP, which should be added to 
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Figure 5. FLI using Cartesian equations. The 1:2 resonance for e = 0.2, 
i = 10°, lo = 0°, = 0°, under the effects of Earth’s oblateness (top left), 

adding Moon + Sun + SRP with A/m = 0.01 [m 2 /kg] (top middle). The 
2:3 resonance for e = 0.3, i = 10°, lo — 0°, Q = 0°, under the effects of 
Earth’s oblateness + Moon + Sun + SRP with A/m = 0.01 [m 2 /kg\ (top 
right). FLI for the 1:3 resonance with e = 0.005, i = 70°, lo — 0°, O = 0°, 
under the effects of Earth’s oblateness (bottom left), adding Moon + Sun 
+ SRP with A/m = 0.01 [m 2 /kg\ (bottom middle), e = 0.3, i = 25°, 
lo — 0°, Q — 0°, under the effects of Earth’s oblateness + Moon + Sun + 
SRP with A/m = 0.01 [m 2 /kg] (bottom right). 


the Hamiltonian (2.2), can be written as (compare with (2.1)) 

Rsrp = ~C r P r Og (—) t - 7 , (6-1) 

m |r — r s | 

where we can use the expansion in terms of the Legendre polynomials: 



with Q being the angle between the Sun and the geocentric radius of the debris (the zero¬ 
th order term T- in the expansion can be neglected, since it does not contribute to the 
Hamiltonian). As for the position of the Sun, we use the following formulae borrowed from 
[18], and normalized with respect to the geostationary distance a geo = 42164.1696 km as 
well as with a unit of time r chosen such that the period of Earth’s rotation becomes 
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equal to 2n: 

r s = 23.71681950544094 (149.619 - 2.499 cos M s - 0.021 cos2 M s ) 

2t r 

M s = -(357.5256 + 35999.04944 t) 

b 360 V ; 

27t 27t 6892 72 

A s = - 282.94 + Ms H-(-sin Ms 4-sin 2 Ms) 

b 360 360 V 3600 3600 ; 

365.242196 # 

36525 • 366.242196 • 2vr ’ 

where Ms is the Sun’s mean anomaly and A s is the ecliptic longitude. Casting together 
these formulae, by the algebraic manipulator Mathematica© we compute the expansion 
of (6.1) up to the third order of the Legendre polynomials, taking care of neglecting in 
the expansion those terms which are multiplied by coefficients less than a specific error, 
whose dimension is kg/m 2 d/r 2 , and that we fix (after trials and errors) equal to 10~ 9 . 
Finally, we average over the mean anomaly to obtain the following expression V///!// for 
the approximate expansion of the potential describing the SRP: 

A / 

Vf/P = a e — ( - 4.838 10” 7 sin(-0.00546061# + u) 
p m V 

- 4.836 10" 7 sin(—0.00546061# + u) - 0.000028751 sin(-0.0027303# + u) 

+ 1.239 10" 6 sin(0.0027303# + u) + 5.425 10“ 6 cos(-0.0027303# + u) 

+ 1.14110~ 7 cos(o; + 0.0027303#)) . 

Clearly, the dimension of is d 2 /r 2 . 

In Figures 6, 7, 8, the FLI is computed as a function of the resonant angle and semi¬ 
major axis, respectively for the 1:2, 1:3, 2:3 resonances. The canonical equations include 
the geopotential and the perturbing effect of the solar radiation pressure. The patterns 
are obtained by using the following area-to-mass ratios: A/m = 0, 1, 5, 15 [m 2 /kg\, and 
are represented in each figure, respectively, in the top left, top right, bottom left and 
bottom right panels. 

By analyzing these results, the main conclusion is the following. Even if the amplitude 
of the resonance is of the order of a few tens of kilometers (top left plots), the solar radi¬ 
ation pressure effect for objects with high area-to-mass ratios yields a web of secondary 
resonances which covers an area of several hundred kilometers. For example, in the case 
of the 1:3 resonance, the amplitude of the resonance for objects with small area-to-mass 
ratio located in planar eccentric orbits with e = 0.3 is about 30 kilometers (Figure 7 top 
left). For A/m = 15 [m 2 /kg] the web of resonances extends on a region of more than 500 
kilometers (Figure 7 bottom right). 
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Figure 6. FLI for the 1:2 resonance, under the effects of the 
geopotential and SRP, for i — 0°, e = 0.1, u — 0°, Q = 0°: 

A/m = 0 [ m 2 /kg\ (top left); A/m = 1 [m 2 /kg\ (top right); 

A/m = 5 [m 2 /kg] (bottom left); A/m = 15 [m 2 /kg] (bottom right). 


7. Two case studies: XMM-Newton and Integral 

As we mentioned in Section 1, two space missions are directly connected with external 
resonances. Precisely, XMM-Newton has a semimajor axis located at the 1:2 resonance, 
while the semimajor axis of Integral is at the 1:3 resonance. However, we must extend 
the analysis performed in the previous sections, since the orbital eccentricities of such 
missions are quite large: XMM-Newton has an eccentricity equal to e = 0.776, while 
it is e = 0.824 in the case of Integral. Such high eccentricities lead to consider some 
different terms in the expansion of the resonant parts associated to these two resonances. 
Moreover, the functions G npq will be developed up to a larger order in eccentricity (14 t/l 
order in the forthcoming computations). In particular, we take N — 3 and among all 
harmonic terms up to degree and order n = m = 3, listed in Table 4, we consider in the 
computations just the most important ones. These terms are reported in bold and their 
explicit expressions are given in Appendix A. 
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Figure 7. FLI for the 1:3 resonance, under the effects of the 
geopotential and SRP, for i = 0°, e = 0.3, u = 0°, Q = 0°: 

A/m = 0 [m 2 /kg\ (top left); A/m = 1 [m 2 /kg\ (top right); 

A/m = 5 [m 2 /kg] (bottom left); A/m — 15 [m?/kg] (bottom right). 


j-t 

N 

terms 

1:2 

4 

72202, 72214, 72226, 7310 1, 73111, 73123, 73201, 73213, 73303, 

72100 , 72112 , 72124 , 73135 , 73225, %237, %315 , 73327, >3339 

1:3 

4 

72204, 72216, 72228, 73100, 73112, 

72101, 72113 , 72125 , 73124 , 73136 , 73203 , 73215 , 73227 , 73239 , 73306 , 73318 , 733210 , 733312 


Table 4. Terms for the 1:2 and 1:3 resonances in the case of large 
eccentricities; the sum of these terms provides R^th' The most relevant 
ones, which are considered in our computations and whose explicit 
expressions are given in Appendix A, are reported in bold. 


The dominant terms for large eccentricities and the amplitudes of the resonances are 
shown in Figure 9 for the 1:2 and 1:3 resonances. For a fast identification, we also 
mark with bullets the positions of XMM-Newton (within the plots referring to the 1:2 
resonance) and Integral (within those referring to the 1:3 resonance). 

Concerning the secular part, since | J30I <C J20 we neglect all terms of order J30. 
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Figure 8. FLI for the 2:3 resonance, under the effects of the 
geopotential and SRP, for i = 0°, e = 0.2, cu = 0°, Q = 0°: 

A/m = 0 [m 2 /kg] (top left); A/m = 1 [m 2 /kg] (top right); 

A/m = 5 [m 2 /kg] (bottom left); A/m = 15 [m 2 /kg] (bottom right). 

Let us first analyze XMM-Newton, whose parameters are e = 0.776 and i = 65.4° ([26]). 
From Figure 9 we see that the leading term is 72214 with an amplitude of the resonance 
around 60-70 km. The orbital data of Integral correspond to e = 0.824 and i = 52.2° 
([13]). In this case 72216 is the dominant term with an amplitude of the resonance of 
about 70 km. These values are in agreement with the FLI plots shown in Figure 10. 

The stable equilibrium point for the 1:2 resonance corresponds to the value of a 12 = 
2cj + A 22 (modulus 180°) with cu = 93° and A22 = 75°; therefore, the equilibrium is located 
at 81° as shown in the upper left panel of Figure 10. Due to the interplay between the 
terms associated to 7 2 2 and those to J 3 1, J 32 , 733, we have a peculiar pattern with a 
double shape surrounding the main island. 

As for the 1:3 resonance, the location of the equilibrium points is influenced by both 
the white and black terms of Figure 9. More precisely, since cu = 302°, A22 = 75° the 
dominant white term gives rise to an equilibrium point at ai 3 = 3cu + A 2 2 (modulus 180°) 
~ 261°, while the equilibrium point corresponding to the black term is at 013 = 2u ; + A 22 
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FIGURE 9. Dominant terms (upper panels) and amplitudes (lower 
panels) of the 1:2 (left panels) and 1:3 (right panels) resonances. The 
black dots represent XMM-Newton on the left panels and Integral in the 
right panels. The colors of the dominant terms are the following. For the 
1:2 resonance: 7-2202 black, 73m blue, 72214 green; for the 1:3 resonance: 
72204 black, 73112 brown, 7-2216 white. The color bar in the lower panels 
provides the amplitude of the resonance in kilometers. 


(modulus 180°) ~ 319°. Because these two terms are comparable in magnitude, the stable 
equilibrium point is in fact located between 261° and 316°, closer to 261°, as shown in 
the upper right panel of Figure 10. 

The middle panels of Figure 10 show the integration of the Cartesian equations for 
8 000 sidereal days, including the effects of the Sun, the Moon and the solar radiation 
pressure. The comparison between the Hamiltonian and Cartesian plots (upper and 
middle panels in Figure 10) is quite different for the 1:2 and 1:3 resonance. Indeed, for 
the 1:3 resonance we have that the Hamiltonian plot provides a fairly exact location of 
the stable equilibrium as well as a good approximation of the separatrix. Of course, the 
dynamics within the resonance island is affected by the influence of Sun, Moon, SRP 
as well as all other terms providing the geopotential at such high eccentricities. For 
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FIGURE 10. Top left: FLI for the 1:2 resonance for e = 0.776, i = 65.4°, 
u> = 93.3°, fl = 55.5°. Top right: FLI for the 1:3 resonance for e = 0.824, 
i = 52.2°, u) = 302°, f! = 103°. Middle panels: FLI in Cartesian formalism 
for the 1:2 (left) and 1:3 resonance (right). Bottom panels: FLI in 
Cartesian formalism for the 1:2 resonance without the effect of the Moon 
(left) and without the effect of Moon, Sun, SRP (right). 


the 1:2 resonance the results are not comparably good (upper and middle left panels in 
Figure 10), since we observe a distortion and displacement of the separatrix. To highlight 
which effect is dominant for the 1:2 resonance, we have integrated the Cartesian equations 
eliminating the effect of the Moon (bottom left panel in Figure 10) and then eliminating 
also the effects of the Sun and SRP (bottom right panel in Figure 10). The results show 
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that the Moon is definitely important in shaping the dynamics of the 1:2 resonance, while 
Sun and SRP have secondary, although still important, effects. 

8. Conclusions and perspectives 

The dynamics of satellites and debris in external resonances is well approximated by 
the Hamiltonian approach based on the study of geopotential including the secular and 
resonant terms. In particular, elementary mathematical techniques provide a fast and 
reliable procedure to reduce the study to the dominant terms (instead of considering 
long series expansions, thus considerably decreasing the computational effort), to locate 
the equilibrium points and to determine their stability, to measure the amplitude of the 
resonance (see Section 4). All these information provide us with a strong background, 
which makes easier to investigate the dynamics of real bodies, when additional effects, 
like those of Sun, Moon and SRP, are considered. 

It is worth to underline that our analysis refers mainly to the long term evolution of 
the semimajor axis, which is influenced directly by the resonant part of the geopotential. 
Therefore, we focused especially on this aspect of the dynamics. However, one should 
mention that the disturbing functions due to Moon, Sun and solar radiation pressure also 
influence the long term behavior of the semimajor axis, though indirectly. More precisely, 
their associated Fourier expansions, with respect to the orbital elements, do not contain 
resonant harmonic terms and, therefore, the canonical equation describing the evolution 
of the semimajor axis (or L) do not change its form, when these additional perturbations 
are considered. However, the long term evolution of the other orbital elements (equiv¬ 
alently, the other Delaunay variables) are directly affected by these perturbations. In 
particular, the eccentricity and inclination vary slowly over time and therefore they indi¬ 
rectly affect the evolution of the semimajor axis, since the canonical equation describing 
the evolution of the semimajor axis involves all orbital elements as parameters. The 
problem is much more complex when shadow effects are considered, especially for high 
area to mass ratio objects. In conclusion, in order to get a clearer picture of how the 
semimajor axis evolves in time, as a result of the interaction of all disturbing forces, one 
needs to quantify the effects of all perturbations as a function of the orbital elements. 

The alternative approach consists in the numerical integration of the Cartesian equa¬ 
tions of motion. In this case it is straightforward to include the effects of Sun, Moon 
and SRP, at the expenses of a longer computational time needed to have equivalent plots 
as in the Hamiltonian case. Moreover, the Cartesian approach, although providing a 
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validation of the Hamiltonian results, does not allow to provide information about the 


location of the equilibria as well as about the shape and width of the resonant region, as 


it was done instead by implementing the Hamiltonian procedure. 

The Hamiltonian approach was also used to investigate the behavior of objects with 
large area-to-mass ratio by suitably expanding the potential describing the effect of SRP 
(see Section 6). The results show the appearance of a web of secondary resonances ranging 
over an area of several hundreds kilometers. In view of these results, it would be inter¬ 
esting to explore the possibility to play with the area-to-mass ratio to stabilize objects 
within specific resonant islands or rather to move them through the web of resonances 
and chaotic layers, which appear when A/m is large. 

In Section 7 we have considered two specific space missions: XMM-Newton, with a 
semimajor axis at the 1:2 resonance, and Integral, with a semimajor axis at the 1:3 
resonance. In both cases the eccentricity is high and therefore we extended our analysis 
by including different terms in the series expansion of the geopotential. Although for 
smaller eccentricities we noticed that the effect of the Moon was more relevant for the 
1:3 resonance, in the case of larger eccentricities we noticed that the lunar attraction 
provokes drastic changes in the shape of both the 1:2 and 1:3 the resonant islands. This 
result suggests that the role of the Moon, as well as that of the Sun - of course, should be 
carefully analyzed (as it will be done in a later work) as a function of distance, eccentricity 
and inclination, in order to have an effective and complete description of the dynamics. 

Appendix A. Leading terms 

We report below the leading terms of the expansion of the secular part given in (3.2). 
For eccentricities larger than 0.5 (Section 7), we considered just the effect of 72 oio : 

>200-2 = T :2022 — 0 , 



We report below the terms of the expansion of the 1:2 resonance, listed in Table 2 or 
written in bold in Table 4. For eccentricities smaller than 0.5, the terms written below 
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have been considered up to second order in the eccentricity. For eccentricities larger 
than 0.5 (Section 7), we considered just the terms for which we wrote the corresponding 
eccentricity functions G npq up to order lA th \ 

fi E R 2 E J2i(3 . 5e 2 ' 


721 nn — 


-{-sini(l + cosi)(l - sin(cri 2 - A 2 i)| , 
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f 3 . 


9e 2 
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|^(1 - cos 2 f)G 2 !4 cos[2(<7i 2 - 2a; - A 22 )]} , 


H E R e J 22 f 3 
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where 

G 202 = ye 2 - 19.167 e 4 + 12.521 e 6 - 3.953 e 8 + 0.703 e 10 - 0.114 e 12 - 0.014 e 14 , 
G 214 = 4.813 e 4 + 0.806 e 6 + 2.809 e 8 + 2.661 e 10 + 2.941 e 12 + 3.144 e 14 , 

G 2 26 = 0.0889 e 6 + 0.044 e 8 + 0.044 e 10 + 0.036 e 12 + 0.03 e 14 , 

G 3i0 -i = -e + 1.25 e 3 - 0.146 e 5 + 0.08 e 7 + 0.042 e 9 + 0.034 e 11 + 0.027 e 13 , 

G '311 = 3 e + 2.75 e 3 + 5.104 e 5 + 7.234 e 7 + 9.697 e 9 + 12.403 e 11 + 15.335 e 13 , 

G 301 = 5 e - 22 e 3 + 25.292 e 5 - 10.889 e 7 + 2.843 e 9 - 0.27 e 11 + 0.151 e 13 , 

G 323 = 1.917 e 3 + 3.708 e 5 + 5.899 e 7 + 8.358 e 9 + 11.065 e 11 + 13.997 e 13 , 

G 313 = 12.833 e 3 - 0.521 e 5 + 9.898 e 7 + 10.008 e 9 + 13.128 e 11 + 16.009 e 13 , 

G 303 = 40.75 e 3 - 161.063 e 5 + 217.8 e 7 - 141.946 e 9 + 55.715 e 11 - 14.141 e 13 , 


and 


a 12 = 2M - 6 + 2u + . 


We report below the terms of the expansion of the 1:3 resonance, listed in Table 2 or 
written in bold in Table 4. For eccentricities smaller than 0.5, the terms written below 
have been considered up to second order in the eccentricity. For eccentricities larger 
than 0.5 (Section 7), we considered just the terms for which we wrote explicitly the 
corresponding eccentricity functions G npq up to order 14 th : 
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Ta lo—i = 
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G 20 4 = 33.3125 e 4 - 86.4188 e 6 + 81.349 e 8 - 40.741 e 10 + 12.781 e 12 - 2.822 e 14 , 

G 216 = 9.897 e 6 - 3.905 e 8 + 5.445 e 10 + 2.201 e 12 + 3.288 e 14 , 

G 228 = 0.163 e 8 + 0.027 e 10 + 0.067 e 12 + 0.048 e 14 , 

G 300 = 1 - 6 e 2 + 6.609 e 4 - 1.953 e 6 + 0.46 e 8 + 0.051 e 10 + 0.07 e 12 + 0.055 e 14 , 

G 3 i 2 = 6.625 e 2 + 2.437 e 4 + 6.876 e 6 + 8.722 e 8 + 11.365 e 10 + 14.181 e 12 + 17.218 e 14 , 


and 


<7 13 = 3 M — 9 + 3u + . 


We report below the terms of the expansion of the 2:3 resonance, listed in Table 2: 
T 2201 = flER / 22 {l(l + cosi) 2 y cos(a- 2 3 -u - 2A 22 )} , 

T 3200 = y sini(l + cosi) 2 (l - 6e 2 ) sin(cr 23 - 2A 32 )| , 

/i E i?|J 32 f 15 . ., o 2 53e 2 . , 1 

73212 = - 2 - i UT 8111 *! 1 - 2 cos* - 3 cos 1 ) —— sm(cr 23 - 2u - 2A 32 ) \ , 

a 4 l 8 8 j 

and 

(T 23 = 3M — 2 9 + 3c 0 T 2fl. 
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